"""
Ideals in Tate algebra

The module gives tools for manipulating ideals in Tate algebras
and, in particular, provides an implementation of the Buchberger
algorithm in this context.

AUTHORS:

- Xavier Caruso, Thibaut Verron (2018-09): Buchberger algorithm

- Xavier Caruso, Thibaut Verron (2019-04): F5-type algorithms (PoTe and VaPoTe)
"""

# ***************************************************************************
#    Copyright (C) 2018 Xavier Caruso <xavier.caruso@normalesup.org>
#                       Thibaut Verron <thibaut.verron@gmail.com>
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 2 of the License, or
#    (at your option) any later version.
#                  https://www.gnu.org/licenses/
# ***************************************************************************

from sage.rings.ideal import Ideal_generic

from sage.structure.richcmp import op_EQ, op_NE, op_LT, op_GT, op_LE, op_GE

from sage.rings.polynomial.polydict cimport PolyDict

from sage.rings.tate_algebra_element cimport TateAlgebraTerm
from sage.rings.tate_algebra_element cimport TateAlgebraElement
from heapq import heappush, heappop

from cysignals.signals cimport sig_check

class TateAlgebraIdeal(Ideal_generic):
    r"""
    Initialize a class for ideals in a Tate series algebra

    EXAMPLES::

        sage: R = Zp(3, prec=10, print_mode="digits")
        sage: A.<x,y> = TateAlgebra(R)
        sage: f = 3*x^2 + 5*x*y^2
        sage: g = 5*x^2*y + 3
        sage: I = A.ideal([f,g]); I
        Ideal (...0000000012*x*y^2 + ...00000000010*x^2, ...0000000012*x^2*y + ...00000000010) of Tate Algebra in x (val >= 0), y (val >= 0) over 3-adic Field with capped relative precision 10

    """

    #@cached_method
    def groebner_basis(self, prec=None, algorithm='VaPoTe', **options):
        r"""
        Compute a Groebner basis of the ideal

        INPUT:

        - ``prec`` -- an integer or ``None`` (default: ``None``), the precision
          at which the computations are carried. If ``None``, defaults to the
          algebra precision cap

        - ``algorithm`` -- a string (default: ``VaPoTe``), the algorithm to
          use in the calculations; available algorithms are:

          - ``buchberger``: classical Buchberger algorithm

          - ``buchberger-integral``: first computes a Groebner basis of the
            ideal generated by the same generators over the ring of integers
            (provides better numerical stability)

          - ``PoTe``: a F5-type algorithm where signatures are ordered by
            position over term

          - ``VaPoTe``: a F5-type algorithm where signatures are ordered
            by valuation over position over term

          We refer to [CVV2019]_ and [CVV2020]_ for a detailed description
          of these algorithms.

        - ``options`` -- extra arguments that are passed in to the
          algorithm; this notably include the keyword ``verbose`` (only
          available for ``PoTe`` and ``VaPoTe``) which is an integer
          defining the verbosity level:

          - ``0``: no verbosity (quiet)

          - ``1``: print each new generator and a notification each time a
            J-pair is popped

          - ``2``: in addition, print the outcome of the treatment of a J-pair

          - ``3``: in addition, print all added J-pairs

          - ``4``: print entire series instead of only their leading terms

        OUTPUT:

        The Groebner basis `(g_1,\dots,g_n)` of this ideal, uniquely determined
        by the following conditions::

        - it is minimal, in the sense that the leading coefficient of `g_i`
          does not divide the leading coefficient of `g_j` if `i \neq j`,

        - it is reduced, in the sense that each term of `g_i` is not divisible
          by leading term of `g_j` for `j \neq i` or the leading term of
          `\pi g_i` where `\pi` is the uniformizer,

        - it is normalized so that the leading coefficient of each `g_i` is
          a power of the uniformizer and moreover, if we are working over a Tate
          algebra (and not its ring of integers), each `g_i` has valuation `0`,

        - it is sorted, in the sense that the leading term of `g_i` is greater
          than the leading of `g_{i+1}` for all `i`.

        .. NOTE::

            The result of this method is cached.

        EXAMPLES::

            sage: R = Zp(3, prec=10, print_mode="digits")
            sage: A.<x,y> = TateAlgebra(R)
            sage: f = 3*x^2 + 5*x*y^2
            sage: g = 5*x^2*y + 3
            sage: I = A.ideal([f,g])
            sage: I.groebner_basis()
            [...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
             ...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
             ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>)]

        The algorithm ``buchberger`` is faster than ``buchberger-integral``
        but may lose more precision::

            sage: R = Zp(2, 5, print_mode='digits')
            sage: A.<x,y> = TateAlgebra(R)
            sage: f = x^2*y^6 + x^4 + 25*y^2 + 2*x^3*y^3 + 10*x*y^4 + 10*x^2*y
            sage: g = x^4*y^5 + x^5*y^2 + x^4 + 5*x^2*y + 2*x^5*y^4 + 2*x^6*y + 6*x^3*y^3
            sage: h = 2*x^6*y^4 + 2*x^4 + 4*x^5*y^2 + 8*x^8*y^2 + 8*x^7*y^3 + 8*x^6*y
            sage: I = A.ideal([f,g,h])
            sage: I.groebner_basis(algorithm="buchberger-integral")
            [...0001*x^4 + O(2^4 * <x, y>),
             ...0001*x^2*y + O(2^4 * <x, y>),
             ...0001*y^2 + O(2^4 * <x, y>)]
            sage: I.groebner_basis(algorithm='buchberger')
            [...001*x^4 + O(2^3 * <x, y>),
             ...001*x^2*y + O(2^3 * <x, y>),
             ...001*y^2 + O(2^3 * <x, y>)]

        TESTS::

            sage: I.groebner_basis(algorithm="F4")
            Traceback (most recent call last):
            ...
            NotImplementedError: available algorithms are 'buchberger', 'buchberger-integral', 'PoTe' and 'VaPoTe'

        """
        if prec is None:
            prec = self.ring().precision_cap()
        if algorithm == "buchberger":
            return groebner_basis_buchberger(self, prec, False, **options)
        elif algorithm == "buchberger-integral":
            return groebner_basis_buchberger(self, prec, True, **options)
        elif algorithm == "pote" or algorithm == "PoTe":
            return groebner_basis_pote(self, prec, **options)
        elif algorithm == "vapote" or algorithm == "VaPoTe":
            return groebner_basis_vapote(self, prec, **options)
        else:
            raise NotImplementedError("available algorithms are 'buchberger', 'buchberger-integral', 'PoTe' and 'VaPoTe'")

    def _contains_(self, x):
        r"""
        Return ``True`` if ``x`` lies in this ideal

        INPUT:

        - ``x`` -- a Tate series

        EXAMPLES::

            sage: R = Zp(3, 10)
            sage: A.<x,y> = TateAlgebra(R)
            sage: f = 3*x^2 + 5*x*y^2
            sage: g = 5*x^2*y + 3
            sage: I = A.ideal([f,g])
            sage: f in I  # indirect doctest
            True
            sage: (f+g) in I  # indirect doctest
            True
            sage: (f+1) in I  # indirect doctest
            False

        TESTS::

            sage: I.random_element() in I
            True

        """
        rgb = self.groebner_basis()
        return (x % rgb).is_zero()

    def _contains_ideal(self, I):
        r"""
        Return ``True`` if ``I`` is contained in this ideal

        INPUT:

        - ``I`` -- an ideal in a Tate series algebra

        EXAMPLES::

            sage: R = Zp(3,prec=10,print_mode="digits")
            sage: A.<x,y> = TateAlgebra(R)
            sage: f = 3*x^2 + 5*x*y^2
            sage: g = 5*x^2*y + 3
            sage: I = A.ideal([f,g])
            sage: A.ideal([f]) < I  # indirect doctest
            True
            sage: I < A.ideal([f])  # indirect doctest
            False
            sage: A.ideal([1]) < I  # indirect doctest
            False
            sage: I < A.ideal([1])  # indirect doctest
            True
        """
        self.groebner_basis()
        return all(f in self for f in I.gens())

    def _richcmp_(self, other, op):
        r"""
        Compare this ideal with ``other`` for the rich comparison
        operator ``op``

        INPUT:

        - ``other`` -- an ideal in a Tate series algebra

        - ``op`` -- a comparison operator

        EXAMPLES::

            sage: R = Zp(3, 10)
            sage: A.<x,y> = TateAlgebra(R)
            sage: f = 3*x^2 + 5*x*y^2
            sage: g = 5*x^2*y + 3
            sage: I = A.ideal([f,g])
            sage: A.ideal([f]) < I
            True
            sage: I < A.ideal([f])
            False
            sage: A.ideal([1]) < I
            False
            sage: I < A.ideal([1])
            True
            sage: I <= A.ideal([f,g])
            True
            sage: I == A.ideal([f,g])
            True
            sage: I <= A.ideal([f])
            False
            sage: A.ideal([f]) <= I
            True
            sage: A.ideal([f]) == I
            False

        """
        if op == op_GT:
            return self._contains_ideal(other) and not other._contains_ideal(self)
        elif op == op_GE:
            return self._contains_ideal(other)
        elif op == op_EQ:
            return self._contains_ideal(other) and other._contains_ideal(self)
        elif op == op_NE:
            return not(self._contains_ideal(other) and other._contains_ideal(self))
        elif op == op_LE:
            return other._contains_ideal(self)
        elif op == op_LT:
            return other._contains_ideal(self) and not self._contains_ideal(other)

    def is_saturated(self):
        r"""
        Return ``True`` if this ideal is saturated.

        The ideal `I` is saturated if `\pi f \in I` implies `f \in I`
        for any `f` in the underlying ring. Here `\pi` denotes a
        uniformizer of the field of coefficients.

        .. NOTE::

            All ideals are saturated when `\pi` is invertible.

        EXAMPLES:

        Over classical Tate algebras (where `\pi` is invertible), this
        method always returns ``True``::

            sage: R = Zp(3, prec=10, print_mode="digits")
            sage: A.<x,y> = TateAlgebra(R)
            sage: f = 3*x^2 + 5*x*y^2
            sage: g = 5*x^2*y + 3
            sage: A.ideal([f,g]).is_saturated()
            True
            sage: A.ideal([f,3*g]).is_saturated()
            True

        The test is only relevant over the rings of integers of Tate
        algebras::

            sage: Ao = A.integer_ring()
            sage: Io = Ao.ideal([f,g])
            sage: Io.is_saturated()
            False
            sage: Io.groebner_basis()
            [...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
             ...0000000001*x*y^2 + ...1210121020*x^2 + O(3^10 * <x, y>),
             ...0000000010*x^3 + ...2222222220*y + O(3^10 * <x, y>),
             ...0000000010*y^2 + ...2101210200*x + O(3^10 * <x, y>)]

        Principal ideals are not always saturated::

            sage: Ao.ideal([3*f]).is_saturated()
            False

        """
        if self.ring().base_ring().is_field():
            return True
        gb = self.groebner_basis()
        for g in gb:
            if g.valuation() > 0:
                return False
        return True

    def saturate(self):
        r"""
        Return the ideal obtained by saturating this ideal.

        In other words, the result is the ideal

        .. MATH::

            (I:\pi^\infty) = \{f \in A : \exists n \in \mathbb{N}, \pi^n f \in I\}`

        where `A` is the underlying ring and `\pi` is the uniformizer of the
        field of coefficients.

        .. NOTE::

            When `\pi` is invertible in `A`, all ideals are saturated.

        EXAMPLES:

        Over classical Tate algebras (where `\pi` is invertible), this
        method always returns the same ideal::

            sage: R = Zp(3, prec=10, print_mode="digits")
            sage: A.<x,y> = TateAlgebra(R)
            sage: f = 3*x^2 + 5*x*y^2
            sage: g = 5*x^2*y + 3
            sage: I = A.ideal([f,g]); I
            Ideal (...0000000012*x*y^2 + ...00000000010*x^2, ...0000000012*x^2*y + ...00000000010)
             of Tate Algebra in x (val >= 0), y (val >= 0) over 3-adic Field with capped relative precision 10
            sage: I.saturate()
            Ideal (...0000000012*x*y^2 + ...00000000010*x^2, ...0000000012*x^2*y + ...00000000010)
             of Tate Algebra in x (val >= 0), y (val >= 0) over 3-adic Field with capped relative precision 10

            sage: I.saturate() == I
            True

        However, the result might be different over the ring of integers
        of a Tate algebra::

            sage: Ao = A.integer_ring()
            sage: Io = Ao.ideal([f,g])
            sage: Ios = Io.saturate(); Ios
            Ideal (...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
                   ...0000000001*x*y^2 + ...1210121020*x^2 + O(3^10 * <x, y>),
                   ...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
                   ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>))
             of Integer ring of the Tate Algebra in x (val >= 0), y (val >= 0) over 3-adic Field with capped relative precision 10

            sage: Io == Ios
            False
            sage: Ios.is_saturated()
            True

        TESTS::

            sage: Io < Ios
            True
            sage: 3*Ios < Io
            True

        """
        if self.ring().base_ring().is_field():
            return self
        gb = self.groebner_basis()
        gens = [ g.monic() for g in gb ]
        return self.ring().ideal(gens)



# Grobner bases computations
############################

# Buchberger algorithm

def groebner_basis_buchberger(I, prec, py_integral):
    r"""
    Compute a Groebner basis of the Tate algebra ideal I using Buchberger's algorithm

    INPUT:

    - ``I`` - an ideal in a Tate series algebra

    - ``prec`` - the related precision at which the initial generators
      are truncated

    - ``integral`` -- a boolean; if ``True``, first compute a
      Grobner basis of the ideal generated by the same generators over
      the ring over the ring of integers

    .. NOTE::

        This function is not meant to be called directly, but through the
        ``groebner_basis`` method of Tate algebra ideals.

    EXAMPLES::

        sage: R = Zp(3, prec=10, print_mode="digits");
        sage: A.<x,y> = TateAlgebra(R)
        sage: f = 3*x^2 + 5*x*y^2
        sage: g = 5*x^2*y + 3
        sage: I = A.ideal([f,g]); I
        Ideal (...0000000012*x*y^2 + ...00000000010*x^2, ...0000000012*x^2*y + ...00000000010) of Tate Algebra in x (val >= 0), y (val >= 0) over 3-adic Field with capped relative precision 10
        sage: I.groebner_basis()  # indirect doctest
        [...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
         ...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
         ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>)]

    """
    cdef list gb, rgb, indices, ts, S = [ ]
    cdef int i, j, l
    cdef TateAlgebraTerm ti, tj, t
    cdef TateAlgebraElement f, g, r, s
    cdef bint do_reduce = True
    cdef bint integral = py_integral

    gb = [ ]
    l = 0
    for f in I.gens():
        if not f:
            continue
        g = f.add_bigoh(f.valuation() + prec)
        if not g:
            continue
        gb.append(g)
        l += 1
    indices = list(range(l))

    # We minimize the family of generators
    rgb = gb[:]
    i = 0
    while i < len(rgb):
        ti = (<TateAlgebraElement>rgb[i])._terms_c()[0]
        for j in range(len(rgb)):
            tj = (<TateAlgebraElement>rgb[j])._terms_c()[0]
            if j != i and tj._divides_c(ti, integral):
                del rgb[i]
                del indices[i]
                break
        else:
            i += 1

    # We compute the initial S-polynomials
    for i in range(l):
        ti = (<TateAlgebraElement>gb[i])._terms_c()[0]
        for j in range(i+1, l):
            tj = (<TateAlgebraElement>gb[j])._terms_c()[0]
            if not ti.is_coprime_with(tj):
                s = (<TateAlgebraElement>gb[i])._Spoly_c(<TateAlgebraElement>gb[j])
                if not s.is_zero():
                    t = s._terms_c()[0]
                    heappush(S, (t._valuation_c(), t._exponent, i, j, s))

    # Main loop of Buchberger algorithm
    # Loop invariant:
    # the S-polynomials of pairs of elements in rgb
    # all reduce to zero modulo (rgb,S)
    while S:
        sig_check()
        # We reduce the Grobner basis if needed
        if do_reduce:
            do_reduce = False
            for i in range(len(rgb)-1, -1, -1):
                g = rgb[i]
                rgb[i] = g._positive_lshift_c(1)
                _, rgb[i] = g._quo_rem_c(rgb, False, True, True)
                gb[indices[i]] = rgb[i]

        # We pop a new S-polynomial
        _, _, i, j, f = heappop(S)
        if i >= 0 and (gb[i] is None or gb[j] is None):
            continue
        _, r = f._quo_rem_c(rgb, False, True, integral)
        if r.is_zero():
            continue

        # We add it to our Grobner basis
        tj = r._terms_c()[0]
        j = len(gb)
        for i in range(j):
            g = gb[i]
            if g is None: continue
            ti = g._terms_c()[0]
            if not ti.is_coprime_with(tj):  # first Buchberger criterium
                s = g._Spoly_c(r)
                if not s.is_zero():
                    t = s._terms_c()[0]
                    heappush(S, (t._valuation_c(), t._exponent, i, j, s))
        gb.append(r)

        # We minimize the Grobner basis
        i = 0
        while i < len(rgb):
            ti = (<TateAlgebraElement>rgb[i])._terms_c()[0]
            if tj._divides_c(ti, integral):
                if indices[i] >= l:
                    heappush(S, (ti._valuation_c(), ti._exponent, -1, -1, rgb[i]))
                    gb[indices[i]] = None
                del rgb[i]
                del indices[i]
            else:
                i += 1
        rgb.append(r)
        indices.append(j)
        # and reduce it
        do_reduce = True

    base = I.ring().base_ring()
    if base.is_field():
        if integral:
            # We need to minimize and reduce the Groebner basis again
            i = 0
            while i < len(rgb):
                ti = (<TateAlgebraElement>rgb[i])._terms_c()[0]
                for j in range(len(rgb)):
                    tj = (<TateAlgebraElement>rgb[j])._terms_c()[0]
                    if j != i and tj._divides_c(ti, False):
                        del rgb[i]
                        break
                else:
                    rgb[i] = rgb[i].monic()
                    i += 1
            for i in range(len(rgb)):
                g = rgb[i]
                rgb[i] = g._positive_lshift_c(1)
                _, rgb[i] = g._quo_rem_c(rgb, False, True, True)
        else:
            rgb = [ g.monic() for g in rgb ]
    else:
        rgb = [ g * base(g.leading_coefficient().unit_part()).inverse_of_unit() for g in rgb ]

    rgb.sort(reverse=True)
    return rgb


# F5 algorithms

cdef Jpair(p1, p2) noexcept:
    r"""
    Return the J-pair of ``p1`` and ``p2``

    INPUT:

    - ``p1`` -- a pair (signature, series)

    - ``p2`` -- a pair (signature, series)

    TESTS::



    """
    cdef TateAlgebraTerm s1, s2
    cdef TateAlgebraElement v1, v2
    cdef TateAlgebraTerm t, t1, t2, su1, su2, sv1, sv2
    s1, v1 = p1  # we assume that s1 is not None
    s2, v2 = p2
    if v1 == 0 or v2 == 0:
        return
    sv1 = v1.leading_term()
    sv2 = v2.leading_term()
    # TODO: Maybe can be made more efficient when we know the elements are "monic"
    t = sv1.lcm(sv2)
    t1 = t // sv1
    t2 = t // sv2
    su1 = t1*s1
    if s2 is None:
        return su1, t1*v1
    su2 = t2*s2
    if su1 > su2:
        return su1, t1*v1
    elif su2 > su1:
        return su2, t2*v2


cdef TateAlgebraElement regular_reduce(sgb, TateAlgebraTerm s, TateAlgebraElement v, stopval) noexcept:
    r"""
    Return the result of the regular reduction of the pair ``(s,v)`` by ``sgb``

    INPUT:

    - ``sgb`` -- a list of pairs (signature, series), candidate to be a strong
      Groebner basis; all leading coefficients are assumed to be a power of the
      uniformizer

    - ``s`` -- a term (the signature)

    - ``v`` -- a series

    - ``stopval`` -- the valuation until which reductions are performed

    TESTS::

        sage: cython(                                                                   # needs sage.misc.cython
        ....: '''
        ....: from sage.rings.tate_algebra_ideal cimport regular_reduce
        ....: def python_regular_reduce(gb, s, v, stopval):
        ....:     return regular_reduce(gb, s, v, stopval)
        ....: ''')

        sage: R = Zp(3, prec=8)
        sage: A.<x,y> = TateAlgebra(R)
        sage: tx = x.leading_term()   # the term x
        sage: ty = y.leading_term()   # the term y
        sage: v = (x + y + 2*x^2*y - x^3).add_bigoh(8)
        sage: p1 = (tx, x^3 + 9*x*y)
        sage: p2 = (ty, x*y + 3*x^2*y)

        sage: python_regular_reduce([p1,p2], tx*ty, v, 8)   # indirect doctest          # needs sage.misc.cython
        (2 + O(3^8))*x^2*y + (1 + O(3^8))*x + (1 + O(3^8))*y + O(3^8 * <x, y>)

        sage: python_regular_reduce([p1,p2], tx, v, 8)      # indirect doctest          # needs sage.misc.cython
        (2 + 2*3 + 2*3^2 + 2*3^3 + 2*3^4 + 2*3^5 + 2*3^6 + 2*3^7 + O(3^8))*x^3
        + (2 + O(3^8))*x^2*y + (1 + O(3^8))*x + (1 + O(3^8))*y + O(3^8 * <x, y>)
    """
    # We assume that the elements of the sgb are such that lt(g) = p^v lm(g) to
    # avoid performing divisions
    cdef dict coeffs = { }
    cdef TateAlgebraElement f
    cdef TateAlgebraTerm lt, factor
    cdef list ltds = [ (<TateAlgebraElement>(d[1]))._terms_c()[0] for d in sgb ]
    cdef list terms = v._terms_c()
    cdef int index = 0
    cdef int i
    cdef bint in_rem

    f = v._new_c()
    f._poly = PolyDict(v._poly.__repn, None)
    f._prec = v._prec
    factor = (<TateAlgebraTerm>terms[0])._new_c()
    while len(terms) > index:
        sig_check()
        lt = terms[index]
        if (not lt) or (lt._valuation_c() >= stopval):
            break
        for i in range(len(sgb)):
            sig_check()
            # The comparison below does not perform a division, it only compares valuation and exponents
            if (<TateAlgebraTerm>ltds[i])._divides_c(lt, integral=True):
                # Here we need the elements of sgb to be "monic"
                factor._exponent = lt._exponent.esub((<TateAlgebraTerm>ltds[i])._exponent)
                factor._coeff = lt._coeff >> (<TateAlgebraTerm>ltds[i])._valuation_c()
                if sgb[i][0] is None or factor * sgb[i][0] < s:
                    f = f - (<TateAlgebraElement>sgb[i][1])._term_mul_c(factor)
                    terms = f._terms_c()
                    index = 0
                    break
        else:
            if lt._exponent in coeffs:
                coeffs[lt._exponent] += lt._coeff
            else:
                coeffs[lt._exponent] = lt._coeff
            del f._poly.__repn[lt._exponent]
            index += 1
    f._poly += PolyDict(coeffs, None)
    f._terms = None
    return f


cdef TateAlgebraElement reduce(gb, TateAlgebraElement v, stopval) noexcept:
    r"""
    Return the result of the reduction of ``v`` by ``gb``

    INPUT:

    - ``gb`` -- a list of reductors

    - ``v`` -- a series

    - ``stopval`` -- the valuation until which reductions are performed

    TESTS::

        sage: cython('''  # optional - sage.misc.cython
        ....: from sage.rings.tate_algebra_ideal cimport reduce
        ....: def python_reduce(gb, v, stopval):
        ....:     return reduce(gb, v, stopval)
        ....: ''')

        sage: R = Zp(3, prec=8)
        sage: A.<x,y> = TateAlgebra(R)
        sage: v = (x + y + 2*x^2*y - x^3*y^2).add_bigoh(8)
        sage: g1 = x*y + 3*x^2*y
        sage: g2 = x^3 + 9*y
        sage: python_reduce([g1,g2], v, 8)   # indirect doctest                         # needs sage.misc.cython
        (1 + O(3^8))*x + (1 + O(3^8))*y + O(3^8 * <x, y>)

        sage: python_reduce([g1,g2], v, 5)   # indirect doctest                         # needs sage.misc.cython
        (1 + O(3^8))*x + (1 + O(3^8))*y + (3^5 + O(3^8))*x^8*y^2
        + (3^5 + 2*3^6 + 2*3^7 + O(3^8))*x^7*y + O(3^8 * <x, y>)
    """
    cdef dict coeffs = { }
    cdef TateAlgebraElement f
    cdef TateAlgebraTerm lt, factor
    cdef list ltds = [ (<TateAlgebraElement>d)._terms_c()[0] for d in gb ]
    cdef list terms = v._terms_c()
    cdef int index = 0
    cdef int i

    f = v._new_c()
    f._poly = PolyDict(v._poly.__repn, None)
    f._prec = v._prec
    while len(terms) > index:
        lt = terms[index]
        if (not lt) or (lt._valuation_c() >= stopval):
            break
        for i in range(len(gb)):
            if (<TateAlgebraTerm>ltds[i])._divides_c(lt, integral=True):
                factor = lt._floordiv_c(<TateAlgebraTerm>ltds[i])
                f = f - (<TateAlgebraElement>gb[i])._term_mul_c(factor)
                terms = f._terms_c()
                index = 0
                break
        else:
            if lt._exponent in coeffs:
                coeffs[lt._exponent] += lt._coeff
            else:
                coeffs[lt._exponent] = lt._coeff
            del f._poly.__repn[lt._exponent]
            index += 1
    f._poly += PolyDict(coeffs, None)
    f._terms = None
    return f


def print_pair(p, verbose):
    r"""
    Return a string representation of the pair ``p``

    INPUT:

    - ``p`` -- a pair (signature, series)

    - ``verbose`` -- an integer

    TESTS::

        sage: from sage.rings.tate_algebra_ideal import print_pair
        sage: R = Zp(3, prec=10, print_mode="digits")
        sage: A.<x,y> = TateAlgebra(R)
        sage: v = 3*x^2 + 5*x*y^2
        sage: s = v.leading_term()
        sage: p = (s,v)

    When ``verbose`` is less than 4, only the leading term of
    the series is printed::

        sage: print_pair(p, 0)
        '(sign = ...0000000012*x*y^2, series = ...0000000012*x*y^2 + ...)'

    When ``verbose`` is at least 4, the entire series is printed::

        sage: print_pair(p, 4)
        '(sign = ...0000000012*x*y^2, series = ...0000000012*x*y^2 + ...00000000010*x^2)'
    """
    if verbose > 3:
        return "(sign = %s, series = %s)" % p
    else:
        s, v = p
        return "(sign = %s, series = %s + ...)" % (s, v.leading_term())

def groebner_basis_pote(I, prec, verbose=0):
    r"""
    Run the PoTe algorithm to compute the Groebner basis of ``I``
    and return it

    INPUT:

    - ``I`` -- an ideal

    - ``prec`` -- the precision at which the Groebner basis is computed

    - ``verbose`` -- an integer (default: ``0``)

    TESTS::

        sage: R = Zp(3, prec=10, print_mode="digits")
        sage: A.<x,y> = TateAlgebra(R)
        sage: f = 3*x^2 + 5*x*y^2
        sage: g = 5*x^2*y + 3
        sage: I = A.ideal([f,g])
        sage: I.groebner_basis(algorithm="PoTe")  # indirect doctest
        [...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
         ...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
         ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>)]

        sage: I.groebner_basis(algorithm="PoTe", verbose=2)  # indirect doctest
        ---
        new generator: ...0000000001*x*y^2 + ...
        0 initial J-pairs
        1 elements in GB before minimization
        1 elements in GB after minimization
        grobner basis reduced
        ---
        new generator: ...0000000001*x^2*y + ...
        1 initial J-pairs
        current J-pair: (sign = ...0000000001*y, series = ...0000000001*x^2*y^2 + ...)
        0 remaining J-pairs
        | new element in SGB: (sign = ...0000000001*y, series = ...0000000010*x^3 + ...)
        | add 2 J-pairs
        current J-pair: (sign = ...000000001*y^2, series = ...0000000010*x^3*y + ...)
        1 remaining J-pairs
        | new element in SGB: (sign = ...000000001*y^2, series = ...0000000010*y^2 + ...)
        | add 3 J-pairs
        current J-pair: (sign = ...000000001*y^3, series = ...0000000010*x^3*y^2 + ...)
        3 remaining J-pairs
        | skip: cover by (sign = ...000000001*y^2, series = ...0000000010*y^2 + ...)
        current J-pair: (sign = ...000000001*x*y^2, series = ...0000000010*x*y^2 + ...)
        2 remaining J-pairs
        | skip: sygyzy criterium; signature = ...0000000001*x*y^2
        current J-pair: (sign = ...000000001*x^2*y^2, series = ...0000000010*x^2*y^2 + ...)
        1 remaining J-pairs
        | skip: sygyzy criterium; signature = ...0000000001*x*y^2
        current J-pair: (sign = ...000000001*x^3*y^2, series = ...0000000010*x^3*y^2 + ...)
        0 remaining J-pairs
        | skip: sygyzy criterium; signature = ...0000000001*x*y^2
        4 elements in GB before minimization
        3 elements in GB after minimization
        grobner basis reduced
        [...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
         ...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
         ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>)]

    We check that :issue:`30101` is fixed::

        sage: I.groebner_basis(algorithm="PoTe", prec=100)  # indirect doctest
        [...0000000001*x^3 + ...2222222222*y + ...000000000*x^2*y^2 + O(3^99 * <x, y>),
         ...0000000001*x^2*y + ...01210121020 + O(3^100 * <x, y>),
         ...0000000001*y^2 + ...01210121020*x + ...0000000000*x^3*y + O(3^99 * <x, y>)]
    """
    cdef TateAlgebraElement g, v
    cdef TateAlgebraTerm s, sv, S, ti, tj
    cdef list terms
    cdef TateAlgebraTerm term_one = I.ring().monoid_of_terms().one()
    cdef bint integral = not I.ring().base_ring().is_field()

    gb = [ ]

    for f in sorted(I.gens()):
        sig_check()
        if not f: # Maybe reduce first?
            continue

        # TODO: this should probably be a single function call
        f = f.monic() << f.valuation()

        if verbose > 0:
            print("---")
            print("new generator: %s + ..." % f.leading_term())
        # Initial strong Grobner basis:
        # we add signatures
        sgb = [ (None, g) for g in gb if g ]
        # We compute initial J-pairs
        l = len(sgb)
        p = (term_one, f.add_bigoh(prec))
        Jpairs = [ ]
        for P in sgb:
            sig_check()
            J = Jpair(p, P)
            if J is not None:
                heappush(Jpairs, J)
        sgb.append(p)

        # For the syzygy criterium
        gb0 = [ g.leading_term() for g in gb ]

        if verbose > 1:
            print("%s initial J-pairs" % len(Jpairs))
        if verbose > 2:
            for s,v in Jpairs:
                print("| " + print_pair((s,v),verbose))

        while Jpairs:
            sig_check()
            s, v = heappop(Jpairs)
            sv = v.leading_term()

            if verbose == 1:
                print("pop one J-pair; %s remaining J-pairs" % len(Jpairs))
            if verbose > 1:
                print("current J-pair: " + print_pair((s,v), verbose))
                print("%s remaining J-pairs" % len(Jpairs))

            # The syzygy criterium
            syzygy = None
            for S in gb0:
                sig_check()
                if S.divides(s):
                    syzygy = S
                    break
            if syzygy is not None:
                if verbose > 1:
                    print("| skip: sygyzy criterium; signature = %s" % syzygy)
                continue

            # We check if (s,v) is covered by
            # the current strong Grobner basis
            cover = None
            for S, V in sgb:
                sig_check()
                if S is not None and S.divides(s):
                    sV = V.leading_term()
                    if (s // S)*sV < sv:
                        cover = (S,V)
                        break
            if cover is not None:
                if verbose > 1:
                    print("| skip: cover by " + print_pair(cover, verbose))
                continue

            # We perform regular top-reduction
            sgb.append((None, v._positive_lshift_c(1)))
            v = regular_reduce(sgb, s, v, prec)
            terms = v._terms_c()
            if terms:
                sv = terms[0]
                if not sv:
                    v = v.add_bigoh(sv._coeff.precision_absolute())
            del sgb[-1]

            if not v or v.valuation() >= prec:
                # We have a new element in (I0:f) whose signature
                # could be useful to strengthen the syzygy criterium
                gb0.append(s)
            else:
                # We update the current strong Grobner basis
                # and the J-pairs accordingly
                vv = v.monic() << v.valuation()
                p = (s,vv)
                if verbose > 1:
                    print("| new element in SGB: " + print_pair(p, verbose))
                count = 0
                for P in sgb:
                    sig_check()
                    J = Jpair(p, P)
                    if J is not None:
                        count += 1
                        if verbose > 2:
                            print("| add J-pair: " + print_pair(J, verbose))
                        heappush(Jpairs, J)
                if verbose > 1:
                    print("| add %s J-pairs" % count)
                sgb.append(p)

        # We forget signatures
        gb = [g for (s,g) in sgb]
        if verbose > 1:
            print("%s elements in GB before minimization" % len(gb))
        if verbose > 3:
            for g in gb:
                print("| %s" % g)

        # We minimize the Grobner basis
        i = 0
        while i < len(gb):
            ti = (<TateAlgebraElement>gb[i])._terms_c()[0]
            for j in range(len(gb)):
                sig_check()
                tj = (<TateAlgebraElement>gb[j])._terms_c()[0]
                if j != i and tj._divides_c(ti, integral):
                    del gb[i]
                    break
            else:
                i += 1
        if verbose > 0:
            if verbose > 1:
                aft = " after minimization"
            else:
                aft = ""
            print("%s elements in GB%s" % (len(gb), aft))
        if verbose > 3:
            for g in gb:
                print("| %s" % g)
        # and reduce it
        for i in range(len(gb)-1, -1, -1):
            sig_check()
            g = gb[i]
            gb[i] = g._positive_lshift_c(1)
            _, gb[i] = g._quo_rem_c(gb, False, True, True)
        if verbose > 1:
            print("grobner basis reduced")
        if verbose > 3:
            for g in gb:
                print("| %s" % g)

    if not integral:
        gb = [ f.monic() for f in gb ]
    gb.sort(reverse=True)
    return gb


def groebner_basis_vapote(I, prec, verbose=0, interrupt_red_with_val=False, interrupt_interred_with_val=False):
    r"""
    Run the VaPoTe algorithm to compute the Groebner basis of ``I``
    and return it

    INPUT:

    - ``I`` -- an ideal

    - ``prec`` -- the precision at which the Groebner basis is computed

    - ``verbose`` -- an integer (default: ``0``)

    - ``interrupt_red_with_val`` -- a boolean (default: ``False``); if
      regular reductions have to be interrupted as soon as the valuation
      raises

    - ``interrupt_interred_with_val`` -- a boolean (default: ``False'');
      if intermediate interreductions of the Groebner basis have to be
      interrupted as soon as the valuation raises

    TESTS::

        sage: R = Zp(3, prec=10, print_mode="digits")
        sage: A.<x,y> = TateAlgebra(R)
        sage: f = 3*x^2 + 5*x*y^2
        sage: g = 5*x^2*y + 3
        sage: I = A.ideal([f,g])
        sage: I.groebner_basis(algorithm="VaPoTe")  # indirect doctest
        [...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
         ...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
         ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>)]

        sage: I.groebner_basis(algorithm="VaPoTe", interrupt_red_with_val=True)  # indirect doctest
        [...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
         ...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
         ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>)]

        sage: I.groebner_basis(algorithm="VaPoTe", verbose=2)  # indirect doctest
        grobner basis reduced
        ---
        new generator: ...0000000012*x*y^2 + ...
        generator reduced
        0 initial J-pairs
        1 elements in GB before minimization
        1 elements in GB after minimization
        grobner basis reduced
        ---
        new generator: ...0000000012*x^2*y + ...
        generator reduced
        1 initial J-pairs
        pop one J-pair; 0 remaining J-pairs
        | add signature for syzygy criterium: ...0000000001*y
        2 elements in GB before minimization
        2 elements in GB after minimization
        grobner basis reduced
        ---
        new generator: ...0000000010*x^3 + ...
        generator reduced
        2 initial J-pairs
        pop one J-pair; 1 remaining J-pairs
        | new element is SGB: (sign = ...000000001*y, series = ...0000000010*y^2 + ...)
        | add 3 J-pairs
        pop one J-pair; 3 remaining J-pairs
        | skip: cover by (sign = ...000000001*y, series = ...0000000010*y^2 + ...)
        pop one J-pair; 2 remaining J-pairs
        | add signature for syzygy criterium: ...000000001*x*y
        pop one J-pair; 1 remaining J-pairs
        | skip: sygyzy criterium; signature = ...0000000001*x^2*y
        pop one J-pair; 0 remaining J-pairs
        | skip: sygyzy criterium; signature = ...0000000001*x^2*y
        4 elements in GB before minimization
        3 elements in GB after minimization
        grobner basis reduced
        [...000000001*x^3 + ...222222222*y + O(3^9 * <x, y>),
         ...0000000001*x^2*y + ...1210121020 + O(3^10 * <x, y>),
         ...000000001*y^2 + ...210121020*x + O(3^9 * <x, y>)]

    We check that :issue:`30101` is fixed::

        sage: I.groebner_basis(algorithm="VaPoTe", prec=100)  # indirect doctest
        [...0000000001*x^3 + ...2222222222*y + ...000000000*x^2*y^2 + O(3^99 * <x, y>),
         ...0000000001*x^2*y + ...01210121020 + O(3^100 * <x, y>),
         ...0000000001*y^2 + ...01210121020*x + ...0000000000*x^3*y + O(3^99 * <x, y>)]
    """
    cdef TateAlgebraElement g, v
    cdef TateAlgebraTerm s, S, sv, ti, tj
    cdef list terms
    cdef bint do_reduce, integral
    term_one = I.ring().monoid_of_terms().one()
    gb = [ ]

    gens = [ ]
    for f in I.gens():
        if f:
            val = f.valuation()
            if val < prec:
                heappush(gens, (val, f.add_bigoh(prec)))

    do_reduce = False
    integral = not I.ring().base_ring().is_field()

    while gens:
        sig_check()
        val, f = heappop(gens)
        if val > prec:
            break

        # We reduce the current Gröbner basis
        if val == 0 or do_reduce:
            for i in range(len(gb)-1, -1, -1):
                sig_check()
                g = gb[i]
                gb[i] = g._positive_lshift_c(1)
                tgtval = val + 1 if interrupt_interred_with_val else prec
                gb[i] = reduce(gb, g, tgtval)
            if verbose > 1:
                print("grobner basis reduced")
            if verbose > 3:
                for g in gb:
                    print("| %s" % g)
            do_reduce = False

        if verbose > 0:
            print("---")
            print("new generator: %s + ..." % f.leading_term())

        tgtval = val + 1 if interrupt_red_with_val else prec
        f = reduce(gb, f, tgtval)
        if verbose > 1:
            print("generator reduced")

        if not f:
            if verbose > 0:
                print("reduction to zero")
            continue

        f = f.monic() << f.valuation()

        if f and f.valuation() > val:
            if verbose > 0:
                print("reduction increases the valuation")
            heappush(gens, (f.valuation(), f))
            continue

        # Initial strong Grobner basis:
        # we add signatures
        sgb = [ (None, g) for g in gb if g ]
        # We compute initial J-pairs
        l = len(sgb)
        p = (term_one, f.add_bigoh(prec))
        Jpairs = [ ]
        for P in sgb:
            J = Jpair(p, P)
            if J is not None:
                heappush(Jpairs, J)
        sgb.append(p)

        # For the syzygy criterium
        gb0 = [ g.leading_term() for g in gb ]

        if verbose > 1:
            print("%s initial J-pairs" % len(Jpairs))
        if verbose > 2:
            for s,v in Jpairs:
                print("| " + print_pair((s,v),verbose))

        while Jpairs:
            sig_check()
            s, v = heappop(Jpairs)
            sv = v.leading_term()

            if verbose == 2:
                print("pop one J-pair; %s remaining J-pairs" % len(Jpairs))
            if verbose > 2:
                print("current J-pair: " + print_pair((s,v), verbose))
                print("%s remaining J-pairs" % len(Jpairs))

            # The syzygy criterium
            syzygy = None
            for S in gb0:
                sig_check()
                if S.divides(s):
                    syzygy = S
                    break
            if syzygy is not None:
                if verbose > 1:
                    print("| skip: sygyzy criterium; signature = %s" % syzygy)
                continue

            # We check if (s,v) is covered by
            # the current strong Grobner basis
            cover = None
            for S, V in sgb:
                sig_check()
                if S is not None and S.divides(s):
                    sV = V.leading_term()
                    # FIXME: should be a monic division
                    if (s // S)*sV < sv:
                        cover = (S,V)
                        break
            if cover is not None:
                if verbose > 1:
                    print("| skip: cover by " + print_pair(cover, verbose))
                continue

            # We perform regular top-reduction
            sgb.append((None, v._positive_lshift_c(1)))
            tgtval = val + 1 if interrupt_red_with_val else prec
            v = regular_reduce(sgb, s, v, tgtval)
            terms = v._terms_c()
            if terms:
                sv = terms[0]
                if not sv:
                    v = v.add_bigoh(sv._coeff.precision_absolute())
            del sgb[-1]

            if v:
                v = v.monic() << v.valuation()
            if v.valuation() > val:
                # We have a new element in (I0:f) whose signature
                # could be useful to strengthen the syzygy criterium
                if verbose > 1:
                    print ("| add signature for syzygy criterium: %s" % s)
                gb0.append(s)
                if v:
                    heappush(gens, (v.valuation(), v))
            else:
                # We update the current strong Grobner basis
                # and the J-pairs accordingly
                p = (s,v)
                if verbose > 1:
                    print("| new element is SGB: " + print_pair(p, verbose))
                count = 0
                for P in sgb:
                    sig_check()
                    J = Jpair(p, P)
                    if J is not None:
                        count += 1
                        if verbose > 2:
                            print("| add J-pair: " + print_pair(J, verbose))
                        heappush(Jpairs, J)
                if verbose > 1:
                    print("| add %s J-pairs" % count)
                sgb.append(p)

        # We forget signatures
        # gb = [ v.monic() for (s,v) in sgb ]
        gb = [ v for (s,v) in sgb ]
        if verbose > 1:
            print("%s elements in GB before minimization" % len(gb))
        if verbose > 3:
            for g in gb:
                print("| %s" % g)
        # we minimize the Grobner basis
        i = 0
        while i < len(gb):
            ti = (<TateAlgebraElement>gb[i])._terms_c()[0]
            for j in range(len(gb)):
                sig_check()
                tj = (<TateAlgebraElement>gb[j])._terms_c()[0]
                if j != i and tj._divides_c(ti, integral):
                    del gb[i]
                    break
            else:
                i += 1
        if verbose > 0:
            if verbose > 1:
                aft = " after minimization"
            else:
                aft = ""
            print("%s elements in GB%s" % (len(gb), aft))
        if verbose > 3:
            for g in gb:
                print("| %s" % g)
        # and reduce it
        do_reduce = True

    # We normalize the final Gröbner basis
    for i in range(len(gb)-1, -1, -1):
        sig_check()
        g = gb[i]
        gb[i] = g._positive_lshift_c(1)
        gb[i] = reduce(gb, g, prec)
    if verbose > 1:
        print("grobner basis reduced")
    if verbose > 3:
        for g in gb:
            print("| %s" % g)
    if not integral:
        gb = [ f.monic() for f in gb ]
    gb.sort(reverse=True)

    return gb
